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ABSTRACT 

First  a  simple  practical  procedure  for  approximating  a  stationary 
Gaussian  process  over  a  finite  interval  by  a  trigonometric  polynomial 
with  predetermined  error  is  described. 

The  approximation  is  then  used  to  calculate  the  distribution  of  the 
maximum,  using  a  novel  Monte  Carlo  method  with  a  control  variable  which 
drastically  reduces  the  variance. 

Finally,  the  outlined  approach  is  compared  to  the  moving-average 
technique  and  shown  to  be  superior  for  continuous-time,  narrow-band 
processes. 


1. 


Introduction 


Many  problems  associated  with  stochastic  processes  in  continuous 
time  cannot  be  solved  exactly  by  analytical  methods.  The  essential 
reason  for  this  lies,  in  the  author's  opinion,  in  the  fact  that  these 
stochastic  processes,  even  over  a  finite  time  interval,  consist  of  a 
non-denumerable  infinity  of  dependent  variables.  It  is  on  account  of 
this  that  ever  since  the  advent  of  the  modern  computer,  researchers  using 
stochastic  processes  have  resorted  to  simulation  to  obtain  the  results 
they  required. 

However,  in  spite  of  the  widespread  use  of  simulation,  there  has  been 
an  extraordinary  dearth  of  theoretical  work  in  this  area. 

While  the  fundamental  principles  of  simulation  of  continuous-time 
stochastic  processes  were  formulated  by  Rice  (1954),  the  standard  exposit¬ 
ion  of  the  simulation  approach  based  on  the  spectrum  has  been  given  by 
Shinozuka  and  Jan  (1972)  and  Yang  (1972,1973). 

Since  then,  there  have  been  swift  developments  in  the  area  of  simulat¬ 
ion  of  Markov-type  and  queueing  stochastic  processes,  which  are  used 
extensively  in  Operations  Research.  See  for  instance  the  special  issue 
of  "Operations  Research",  Vol .  31,  No.  6  (Nov.  1983),  dedicated  to  simulat- 


On  the  other  hand  there  have  been  only  a  few  scattered  papers  on 
simulation  of  continuous-time  stochastic-processes  based  on  the  spectrum, 
the  most  easily  accessible  of  which  are  Bi’ly  and  Bukoveczky  (1976), 

Mihailov  (1978)  and  Malyshev  and  Palagin  (1981).  See  also  Bily  and  Cacko 
(1982)  and  Kropa'c  (1981). 

All  proposed  simulation  schemes  rely  on  approximating  the  stochastic 
process  by  a  function  of  a  finite  number  n  of  random  variables.  This 
is  almost  invariably  justified  by  showing  that  the  approximation  converges 
in  some  sense  to  the  original  stochastic  process  as  n  tends  to  infinity. 

However,  this  type  of  justification  gives  no  indication  of  the  size 
of  the  error  committed  when  n  is  taken  to  be  a  comparatively  small  number. 
In  many  applications  the  size  of  this  error  is  not  very  crucial,  because 
the  original  process  is  not  known  with  a  high  precision  anyway. 

But  in  the  application  addressed  to  in  this  paper,  namely  the  numerical 
calculation  of  a  distribution  function,  given  an  exact  formula  for  the 
spectrum,  it  is  of  fundamental  importance  to  know  the  size  of  the  error. 

In  this  paper,  two  new  ideas  are  put  forward  to  obtain  a  computational 
technique  which  is  more  efficient  than  those  previously  proposed,  and  for 
which  exact  bounds  for  the  error  can  be  evaluated. 

The  first  idea  can  be  applied  to  any  stationary  Gaussian  process  X(t) 
with  given  spectral  density  f(x)  and  consists  in  representing  the  process 
over  a  finite  interval  by  a  finite  trigonometric  polynomial  in  such  a  way 
that  the  error  committed  can  be  exactly  bounded. 
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The  second  idea  relates  to  the  calculation  of  the  distribution 
function  of  the  maximum.  The  finite  trigonometric  approximation  is  used 
to  simulate  the  process.  The  maximum  functional  is  then  split  up  into 
two  independent  components.  The  distribution  of  one  component  is  cal¬ 
culated  exactly,  while  that  of  the  other  is  obtained  from  the  simulation. 
The  two  components  are  then  combined,  yielding  a  high  precision  estimate 
for  the  distribution  of  the  maximum. 


The  germ  of  the  last  idea  is  to  be  found  in  Deak  (1980),  but  its 
application  to  the  calculation  of  the  distribution  of  the  maximum  of  a 
stochastic  process  was  given  independently  by  the  author  in  Hasofer  (1982) 
See  also  Moran  (1984). 


The  importance  of  using  well-defined  bounds  for  the  approximation 
error  is  further  highlighted  by  the  results  given  in  Lyon  (1970),  where  it 
is  shown  that  approximations  which  might  be  satisfactory  for  low  levels 
might  be  completely  useless  for  studying  the  behaviour  of  the  stochastic 
process  at  high  levels. 


2.  Construction  of  approximating  polynomial 


In  Hasofer  (1982),  it  is  shown  that,  over  a  finite  interval  (-T,  +T), 
a  real,  zero-mean,  Gaussian  stationary  stochastic  process  X ( t )  with 
spectral  density  f(x)  can  be  approximated  by  a  finite  random  trigonomet¬ 
ric  sum  W^(t),  in  the  sense  that  there  exist  two  stochastic  processes 
W^(t)  and  Y(t),  such  that  X(t)  +  Y(t)  and  WN(t)  +  WN(t)  are 


S-frWS&w- v^v-W-Sv: 


vs1- 


(b)  Choose  a  minimal  subset  of  integers,  N,  such  that 


where  N  is  the  complement  of  N  in  the  set  of  non-negative 
integers. 


▼ 
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(c)  Then 


wN(t) 


wNU) 


Y(t) 


l  a  (yn  COS  ^  +  w  sin 
neN  n  n  2T  n  2T 


.  nut, 


I-  <Vn  TT  +  "n  *’»  TT>  ' 
ncN 

IWn™  If  +  “nsinTf)  ■ 


.  nut. 


where  the  (Vn)  and  the  {W^}  are  sequences  of  independent 
Gaussian  random  variables  with  zero  mean  and  unit  variance. 


The  notion  of  "smallness"  in  this  context  refers  to  the  two  coefficients 
oQ  and  o2  defined  for  the  process  Z(t)  by 


a2o  =  Var  LZ(t) ] 


and 


j2  =  Var  tZ'(t)l 


2 ,  2- 


;ince  Pr{max  |Z(t)  j >x}<2{l-$(x/o  )}  +  (2To~/vo  )exp[-^(x  /0~~\ 
-T<tsT  °  6  ° 


...(2) 


where  ${x)  is  the  standard  normal  cumulative  distribution  function  and 
the  right-hand  side  of  equation  (2)  is  an  increasing  function  of  both  c(. 
and  02»  for  x  >  o2  (see  Hasofer  1982). 


'-\V.  /. V.'--  .V .V  v,'.' 
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It  is  to  be  noted  that  by  taking  N  sufficiently  large,  W^(t) 
can  be  made  as  small  as  required,  while  this  is  not  the  case  with  Y(t). 

The  process  Y(t)  is  derived  from  the  negative  coefficients  of 
the  Fourier  expansion  of  R(t)  over  the  finite  interval  (-2T,  +2T). 

It  is  well  known  that,  as  T  -*•  ®,  the  Fourier  series  tends  to  the  corres¬ 
ponding  Fourier  integral  with  non-negative  integrand.  Thus,  Y(t)  can 
be  made  as  small  as  required,  by  taking  T  sufficiently  large. 

As  an  example  of  the  application  of  the  above  idea,  one  can  approximate 
the  distribution  of  the  maximum  of  X(t)  over  (-T,  +T)  by  that  of 
WN(t)  as  follows: 

Let  Mx  =  max  X(t) 

M1  =  max  ( t ) 

U(t)  =  Y(t)  +  W^(t) 

M  =  max  | U ( t ) | 

where  the  maximum  is  taken  over  all  values  of  t  in  (-T,  +T). 

It  is  then  easy  to  establish  the  following  inequality,  which  holds  for 
every  k. 

P(M1>x+2k)  -2P(M>k)  <  P(Mx>x)  <  P(M]>x-2k)  +2P(M>k) 


(see  Hasofer  (1982)). 
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Clearly  the  error  depends  on  P(.M  >  k),  which  in  turn  is  bounded  as 
in  (2),  Thus  a  measure  of  the  error  is  the  magnitude  of  the  two  coefficients 


Var  CU(t)]  =  [b*  +  La* 
n  n  ncN  n 


3.  Relation  to  discrete  approximation  of  distributions 

Looking  at  the  spectrum  of  the  process  X(t),  which  has  been  assumed  to 
be  continuous  (equation  (1)),  and  that  of  the  process  W^(t),  which  is 
discrete,  when  ( t )  is  continued  to  (-<*>,  +°°)  by  periodicity,  we  see  that 

the  approximation  developed  will  produce  (apart  from  a  scale  factor)  a  dis¬ 
crete  approximation  to  the  absolutely  continuous  probability  density  function 
f(x)/R(0). 

Such  approximations  have  been  shown  to  be  very  useful  in  many  practical 
problems,  for  example  in  the  evaluation  of  moments  of  functions  of  random 
variables.  (See  e.g.  the  much  quoted  paper  of  Rosenlueth  (1975)).  The  proced¬ 
ure  described  in  principle  in  Section  2  and  detailed  in  the  sequel  provides 
a  simple  algorithm  for  obtaining  a  discrete  approximation  with  many  optimal 
properties  which  can  be  derived  from  the  optimality  properties  of  Fourier 
polynomials. 


Mali 


4. 


Optimal  property  of  the  proposed  representation 


Let  W*(t)  be  another  approximation  to  X(.t),  with  covariance 
function  R*(t).  Expanding  R*(t)  in  a  Fourier  Series  over  the  interval 
(-2T,  +21),  we  have,  say 


R*lt)  =  y 3  cos 


nut 
n  w  2T 


n=o 

For  simplicity  let  us  write  en  for  cos(niTt/2T) ,  x+  for  max(0,x)  and 
x"  for  -min(0,x). 


Then 


R(t)  -  I„*en 

n=o 


)a~  e  . 
n 

n=0 


and 


R*(t)  =  ye+  e„  -  y  B"*  e  . 

n=0n  n  n=0  n  n 


It  follows  that 


or 


R(t)  +  Ken  ‘ 

=  R*tt)  ♦  K«;  -sn)+en  - 

R(t)  *  IV„  *  ‘  3  nr%  =  R‘(t)  *  I(“n  '8n,+en- 


This  leads  to  the  representation 

X(t)  +  Y(t)  +  Y*  ( t)  i  W*(t)  +  Z*(t), 

where  Y*(t)  has  covariance  functicn  y(a*  -  3  )  ep  and  Z*(t)  has 

variance  function  y(a+  -  3  )+e  • 

L  n  n  n 


co- 


‘i 

I 


N 

:<i 

y 

y 


.■> 

•1 


I 

s 

i 
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Measuring  the  discrepancy  between  X(t)  and  W*(t)  as  in  the 


previous  section,  we  see  that  both  coefficients 

Var[Y(t)  +  Y*(t)  +  Z*(t)] 

Var[Y'(t)  +  Y*'(t)  +  Z* • (t )  ] 

are  minimized  by  choosing  Sn  =  a*  ,  which  implies  that  W*(t)  has  the 
same  distribution  as  W^(t)  +  W^(t).  Finally  a  standard  least  squares 
argument  shows  that  the  trigonometric  polynomial  of  order  N  which  gives 
the  smallest  error  as  defined  above  is  indeed  WN(t). 

5.  Some  practical  formulae  for  evaluation  of 
the  coefficients  and  the  remainder  terms 


For  the  purpose  of  determining  the  appropriate  set  N,  as  well  as  the 
value  of  T,  it  is  best  to  start  with  the  spectral  density  function  f ( A ) 
of  X(t).  This  is  an  even  function  of  x  and 


RU) 


/■  +  00 

—  CD 


cos(xt)f(x)  dx. 


The  Fourier  coefficients  of  R(t)  in  (-2T,  +2T)  are  given  by 

,2T 


1 


2T 


R  ( t )dt  , 


1 

T  j 


T> 

2T 


R(t)  cos  dt  (n  >  o) , 


and  these  can  be  rewritten  as 


f+”  sin  2xT  ... 


-10- 


a 


n 


sin(2xT-rnr)  sin(2AT+nn ) 

(2xT-rvn)  (2xT+nn) 


f ( x )d> ,  (n  >  o) . 


This  last  formula  can  be  rewritten  more  conveniently  as 


,-fco 

an  *  H>" 

J  —  oo 


4xT  sin  2xT 

2~2 — 2~2  f(x)dx  ,  (n  >  o)  . 


-  (4x  T  -n  tt  ) 


(3) 


One  of  the  most  interesting  cases  is  that  of  a  "narrow  band"  stochastic 
processes,  where  f(x)  is  of  the  form 


f(x) 


_L 

2o 


X  -X 

h(-^-) 

a 


+ 


X  +X 


9 


and  h(x)  is  a  density  function  with  mass  concentrated  in  a  neighbourhood 
of  the  origin. 


Formula  (3)  can  then  be  rewritten  as 


a 

n 


4xT  sin  2xT 
(4x  T  -n  it  ) 


h(z)dz,  n 


o 


with  X  =  oz  +  X  , 

o 


and 


<■+00 


sin  2xT 
2xT 


h(z)dz, 


with  A  as  above. 


1 


I 

-0 


y 
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To  choose  appropriate  values  for  N  and  T,  we  notice  that  the 
distance  between  the  angular  velocities  of  two  successive  terms  such  as 
an  and  an+.j  is  n/2T.  Suppose  that  most  of  the  mass  of  the  spectral 
density  function  on  the  positive  axis  is  concentrated  in  the  interval 
(A0-ko,Ao+ko) ,  where  k  is  some  real  positive  number.  We  now  attempt 
to  approximate  f(x)  by  means  of  2  M+l  masses  equidistantly  placed  on 
the  interval  (xo-ko,xo+ka) .  We  must  then  have 

2ko  _  _jt_ 

2M  2T 

or  2T  =  ■— • 

ko 

The  middle  index  Nq  will  then  be  given  by 


where  [  ]  denotes  as  usual  the  integral  part,  and  the  set  of  indices  chosen 

for  will  be 

N  -M,  N  -M+l,  ....  Nq+M. 


It  will  then  be  necessary  to  check  that  the  values  of  ap  obtained  are  not 
negative.  Any  negative  coefficients  will  have  to  be  removed.  The  values 
of  oQ  and  0£  for  any  of  the  three  processes  W^(t),  W^(t),  Y(t)  will  be 
given  by 


-12- 


where  the  sunmation  extends  over  the  appropriate  set  of  indices. 

In  practice,  provided  fCO  decays  sufficiently  fast  as  \  -*•  «, 

only  a  comparatively  small  number  of  coefficients  will  be  necessary  to 

2  2 

adequately  represent  the  covariance  function,  and  then  aQ  and 

can  be  evaluated  for  each  of  the  three  processes  by  using  only  the  small 

set  of  non-negligible  coefficients. 


As  an  example  of  the  application  of  the  above  methods,  we  consider  a 
narrow-band  process  with  the  spectral  density  function 


A  -A  A  +A 

.  *<-V>  +  ♦(-£->  . 


where  4>  is  the  density  function  of  the  standard  normal  distribution 


We  choose  the  following  parameters 


15tt  , 


2ti 


We  then  choose  2T  =  l  and  k  =  5,  which  gives  us  M  =  2Tk0/Ti  =  10.  The 


index  corresponding  to  xQ  will  then  be  Nq  =  2TAQ/Tt  =  15. 


*.  /.  \ 


3- 


Applying  the  above  fonnulae,  with  numerical  integration  performed 
by  Simpson's  rule  over  a  z-interval  of  (-7,  +7)  with  500  intervals, 
we  obtain  the  coefficients  for  W^(t)  given  in  Table  1. 

Coefficients  beyond  n  =  35  were  negligible,  and  there  were  no 

significant  negative  coefficients,  so  that  Y ( t )  in  this  case  was  neg- 

2  2  2 

ligible.  The  values  of  oQ  and  for  X ( t )  were  oQ  =  1  and 

2  2?  ? 

°2  =  o  +  0  '  222it  =  2260’ 

while  the  corresponding  values  for  W^(t)  were 
aZQ  =  8.58  x  10'7 
and  o g  =  4.98  x  10'3, 

showing  that  the  approximation  obtained  with  twenty  coefficients  is  high 


TABLE  1 


Fourier  Coefficients  for  W^(t) 


a 

n  n 


a  =  a2 
n  n 


5 

7.43 

X 

,o-7 

15 

0.199 

6 

7.99 

X 

to-6 

16 

0.176 

7 

6.69 

X 

to-5 

17 

0.121 

8 

4.36 

X 

to-4 

18 

6.48 

X 

10 

9 

2.22 

X 

to-3 

19 

2.69 

X 

10 

10 

8.76 

X 

to-3 

20 

8.76 

X 

10 

11 

2.70 

X 

CM 

1 

O 

21 

2.22 

X 

10 

12 

6.48 

X 

CM 

1 

O 

22 

4.36 

X 

10 

13 

0.121 

23 

6.69 

X 

10 

14 

0.176 

24 

7.99 

X 

10 
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7.  Calculation  of  P(M^  >  x) 


It  was  pointed  out  in  Hasofer  (1982)  that  the  distribution  of  the 
maximum  of  ( t X  can  be  calculated  directly  by  means  of  a  multiple 

integral.  If  there  are  k  elements  in  N,  the  multiple  integral  extends 
over  the  2k-dimensional  ball  of  unit  radius. 


To  perform  the  integration,  however,  one  cannot  use  efficiently  any 
deterministic  method.  As  was  pointed  out  by  Davis  and  Rabinowitz  (1975), 
the  Monte  Carlo  method  is  the  preferred  approach  for  integration  in  spaces 
of  large  dimensions. 

In  the  light  of  the  preceding  paragraph,  one  can  recast  the  problem  in 
a  purely  probabilistic  framework,  as  follows: 


We  have 


Vt)=  I  anlVn  coslif  *  “n^TT1 

neN 


Write  R2  =  Y  (V2  +  W2).  Then  R2  has  clearly  a  chi-squared  dist- 
neN  n  n 
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V  I  I 


of  the  trigonometric  polynomial 


Find  the  maximum 


TNi>(t)  SlnEff> 


R  ncN 


Estimate  P(M^  >  x)  by 

>  *>  *  5  l  P(R  *  7TTJ 


«;<*>) 


i=l 


Properties  of  the  estimator  P 


Clearly  P  is  an  unbiased  estimator  of  P CM-j  >  x)  since 

E(Pl  ■  Q  J,  E[P(R ’ I  <(1)>]  ■  E  [P<R  >  ^  |M'’] 


A 

For  large  n  ,  P  ,  tends  to  be  normally  distributed,  on  account 
of  the  central  limit  theorem. 


The  variance  of  P  can  easily  be  estimated.  Let 

pi  '  PO  >-*xt7K(,)) 

M1 

Then  an  unbiased  estimator  of  the  variance  of  P.  is 

s  ‘pi>  *  TFTT  J,(pi  '  p)  • 

A 

and  the  corresponding  unbiased  estimator  of  the  variance  of  P  is 


-18- 


By  taking  a  sufficiently  large  number  of  simulations,  the  variance 
of  the  estimator  can  be  made  as  small  as  required. 

2 

The  technique  of  separating  out  R  and  using  its  exact  distribution 
amounts  to  a  "control -variable"  variance  reduction  technique  in  Monte- 
Carlo  integration.  See  e.g.  Rubinstein  (1981). 

9.  Some  numerical  results 


The  technique  described  in  the  preceding  section  was  applied  to  the 

trigonometric  polynomial  given  in  Section  5,  and  the  following  results 
were  obtained. 


(1)  M*: 

Mean:  0.327 


Standard 

deviation: 

0.077 

of  M1 

X 

A, 

P(Mn 

>  x) 

Vs2(pi> 

4 

1.82 

x  10'3 

00 

X 

o 

1 

ro 

5 

9.68 

x  10'6 

2.03  x  10‘4 

6 

2.917 

-9 

x  10 

6.22  x  IQ"8 

(3) 


Number  of  simulations  required. 


Suppose  we  require  a  coefficient  of  variation  of  5%  for  P. 
the  number  of  simulations  required  turns  out  to  be: 


x 

4 


40,000 


175,900 


.  ....... 


182.700 


Then 


v  'i,c 
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These  numbers  are  well  within  the  capabilities  of  an  average  computer. 

(4)  Computer  time  per  simulation 

The  simulation  was  carried  out  both  on  a  large  main  frame  computer  (Control 
Data  Cyber)  and  on  a  small  personal  computer  (IBM  PC  XT).  The  average  time  per 
simulation  for  the  main  frame  was  1.2  seconds,  while  the  average  time  for 
the  IBM  PC  was  18  seconds. 

However,  there  was  no  difficulty  in  obtaining  a  large  number  of  simulations 
on  the  PC,  since  it  could  be  left  running  over  week-ends,  and  for  each  week-end 
about  10,000  simulations  could  be  accumulated. 


Advantage  of  the  above  technique 


One  might  think  that  with  the  above  large  number  of  simulations,  it  might 
not  be  necessary  to  go  through  all  the  steps  of  Section  8.  One  might  be  tempted 
to  simulate  W^(t)  directly  and  estimate  the  tail  probabilities  of  M-j  directly 
from  its  histogram.  To  appreciate  the  large  gain  effected  by  the  procedure  of 
Section  8,  it  is  only  necessary  to  calculate  the  number  of  simulations  necessary 
to  achieve  the  same  coefficient  of  variation  as  in  Section  3,  namely  5»,  by 
direct  simulation. 
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11.  Approximation  to  the  distribution 
of  the  maximum  of  X(t) 

As  mentioned  above  the  following  inequality  holds 

PCM1>x+2k)-2P(M>k)<P(Hx>x)<P(M1>x-2k)+2PCM>k). 

Thus  we  can  approximate  the  distribution  of  the  tail  of  the  maximum  of 
X(t)  over  (-T,  +T)  by  the  distribution  of  the  tail  of  the  maximum  of 
W^(t),  calculated  in  Section  6. 

The  term  P(M  >  k)  can  be  bounded  by  using  formula  (2). 

Suppose  we  take  k  =  0.005.  It  then  turns  out  that 

P(M  >  k)  <  10'5. 

Choose  x  =  3.  It  then  turns  out  that 

P (M 1  >  3.01)  =  0.0624 
P(M]  >  2.99)  =  0.0653 

Thus 

0.0624  ^  P(MX  >  3)  s  0.0654 


corresponding  to  a  margin  of  error  of  ±  3%. 
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12.  The  modulus  of  variation  of  T^(t) 


One  of  the  great  advantages  of  the  trigonometric  polynomial  representat¬ 
ion  is  that  upper  bounds  for  the  variability  of  the  approximation  can  be 

★ 

easly  obtained.  The  approach  will  be  illustrated  by  the  calculation  of 
based  on  the  maximum  of  a  finite  set  of  values  of  t. 


^1  =  max  T^ ( t . ) 

where  the  t^  form  a  mesh  of  span  A  in  (-T,  +  T). 

★ 

Let  tg  be  the  value  of  t  for  which  T^(t)  reaches  ,  and  let 

s  s  Vl  • 

Now  at  tQ  we  must  have  dT^|dt  =0  with  probability  one.  However 

(t  _  t  )2 

W  =  W  +  ^i  "  WV  +  2  v  TN  ’ 


where  t^  <  t  <  tQ. 


W  -  W  - 


(ts  -  tnr  , ,  , 


and  since  obviously  T^(t^)  s  <  T^Uq} 


we  have 


*  **  o 

|M,  -  Mi  I  *  £1  M 


where 


M0  =  max)  -  Tn  (t)| 
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But  -TN  (t)  *  %>2  LnVAncos!lT  *  V1"^  • 

ncN 


so  that 


V?  i'T>'‘«5r>2  I"W»2  ’ 

ncN 


<-  (|f)2(  I  "4#1/2 

LX  ncN  n 


2  2 

using  Schwarz's  inequality  and  the  fact  that  ]_  (A  +  B^)  =  1. 

neN 

In  this  way  we  can  choose  the  span  A  in  advance  so  as  to  make  the 

★ 

error  in  calculating  as  small  as  required. 


13.  Comparison  with  the  moving  average  approach 

The  most  popular  simulation  technique  used  in  recent  times  has  been  the 
moving  average  approach  (see  e.g.  Journel  (197^,  Bily  and  Bukoveczky  (1976) 
and  Journel  and  Huijbregts  (1978)J.  A  brief  (and  somewhat  improved)  presentat- 
i on  foil ows . 

The  method  simulates  a  real,  zero-mean  Gaussian  stationary  stochastic 

process  X(t)  with  given  covariance  function  R(t)  at  a  discrete  sample  of 

Doints  {t  :  nc(-°°,  +°°}  distant  1/2W  from  each  other.  Let  X ( t  )  =  X  and 
r  n  n  n 

R(t  )  =  r  .  The  simulation  attempts  to  represent  a  as  a  moving  average  of 

n  n  n 

independent  standard  normal  variables  (V  '•  nr(-«,  +“>)}  ,  in  the  form 


)  c  v 

L  r  n-  r 

p  —  _  oo 


(4) 


where  the  Cn  are  coefficients  to  be  determined.  It  turns  out  to  be 
convenient  to  set 


C  =  C  . 
-n  n 


r  =  E  { X .  X  }  =  V  C  C  . 
n  k  k+n  L  r  r+n 


Let  now  f^w)  =  \  r  exp(icjn)  (which  is  necessarily  £  0)  be  the 

n=-°° 

power  spectrum  of  Xn, 


and  let 


C(lu)  =  \  C  exp(iwn). 

n=-°°  n 


fn(u>)  =  C(u)CC-w). 


However,  from  the  assumption  C  =  C  ,  we  conclude  that 

-n  n 


C (lu)  =  C(-w)  =  Cq  +  2  y  Cncoston. 

n=l 


It  follows  that  C(w)  can  be  taken  as  /fp(oj),  and  we  then  have,  from  well 
known  Fourier  results 

7T 

C  =  -  I  /fn(w)  cos  on  du>,  n  =  0,1,2,..,  .  (5) 

n  71  j0  U 

With  the  Cp  calculated  from  (5),  Xp  given  by  (4)  has  a  covariance 
function  exactly  equal  to  that  of  X ( t )  at  all  points  tp. 


Note  that  the  variance  of  X„  is  given  by  R ( 0 )  =  )  =  cil  +  2  )  C 

n  _ '■  n  u  _L'i 

n--a-  n- 1 


*1* 
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To  obtain  an  approximation  to  X^  involving  only  a  finite  number  of 
random  variables  we  choose  a  subset  of  integers  N  such  that 


I-  c; 


<  < 


I  c; 


neN  ncN 

where  N  is  the  complement  of  N  in  the  set  of  non-negative  integers. 


We  can  then  write 


where 


and 


X  =  w  +  w 
n  n  n 


W  =  Y  C  X 
n  r  n-r 

reN 


W  =  V  C  X  . 
n  r  n-r 

reN 


For  fixed  n. 


W  and  W  are  independent.  However  for  n  f  m  W  and 
n  n  1  n 

I 

Wm  are  correlated,  unlike  the  decomposition  of  trigonometric  polynomials. 

This  approach  is  in  a  certain  sense  the  dual  of  the  approach  advocated  in  this 
paper,  in  that  the  trigonometric  polynomial  represents  the  original  process 
over  a  finite  time  interval  (-T,  +T),  resulting  in  the  spectrum  being  dis¬ 
crete,  while  in  the  moving  average  process  the  time  parameter  is  discrete, 
resulting  in  a  spectrum  which  is  band-limited. 

One  advantage  of  the  moving  average  process  is  that  the  simulation  can  be 
continued  for  as  long  an  interval  of  time  as  required,  making  it  particularly 
suitable  for  real-time  simulation. 


On  the  other  hand,  the  trigonometric  polynomial  offers  the  advantage  of 
an  analytical  expression  in  closed  form,  which  makes  it  eminently  suitable 
for  applications  where  e.g.  differentiation  and  integration  of  the  process 
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must  be  performed.  An  example  of  the  power  of  this  approach  has  been  given 
in  Section  12.  No  similar  deterministic  error  estimation  can  be  performed 
with  the  moving  average  process. 

But  the  main  advantage  of  the  trigonometric  polynomial  approach  lies  in 
its  ability  to  deal  with  narrow  band  processes.  A  perusal  of  Section  5  shows 
that  for  fixed  T,  the  number  of  terms  of  the  polynomial  decreases  when  the 
band  becomes  narrower.  On  the  other  hand,  when  the  band  becomes  narrower,  the 
covariance  function  decays  more  slowly,  and  more  terms  are  needed  in  the  moving 
average  approach. 

This  property  of  the  moving  average  can  be  illustrated  by  calculating  the 

1 

CnS  for  the  example  given  in  Section  5.  For  points  distant  2W,  we  find 
that  approximately 


12  2 

Cq  -  Cpfexp (--^o^n  )}cosnpQ  > 


where  Oq  =  /2o/2W  and  Pq  =  Xq/2W. 


It  is  clear  that  the  smaller  o,  the  larger  the  number  of  required  coefficients 


For  example,  it  we  take  2W  =  200  (a  number  which  was  actually  used  in 
Section  9)  and  o  =  2tt,  it  turns  out  that  we  need  at  least  80  coefficients 
(as  compared  with  20  coefficients  for  the  trigonometric  polynomial),  and,  of 
course,  for  each  simulation  of  X(t)  for  200  points,  we  need  200  +  (80  x  2) 
=  360  simulations  of  independent  standard  normal  variables,  as  against  only 
40  for  the  trigonometric  polynomial.  It  is  true  that  for  the  trigonometric 
polynomial  we  need  to  calculate  a  large  number  of  sine  and  cosine  values,  but 
these  can  be  efficiently  calculated  by  using  Fast  Fourier  Transform  techniques. 


A  list  of  approximate  values  of  Cn  is  given  in  Table  2  for  the  values 


given  in  the  above  Daraaraoh. 


r  h.  A  *.  *  *  *  ■  *  »  >  v  -  A  -  A  «  *■  A  .  *  -. 
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TABLE  2 


Moving  Average  Coefficients  for  X(t) 
CQ  =  0.2238 


n 

cn 

n 

c 

m 

10 

-0.1434 

60 

2.72  x  10"8 

20 

-9.856  x  10‘9 

70 

-0.00122 

30 

0.0651 

80 

3.53  x  10"4 

40 

-0.0461 

90 

-8.30  x  10'5 

50 

0.0135 

100 

-3.70  x  10"8 
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